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Abstract — The  deformation  of  micro-twinned  polycrystalline  magnesium  (Mg)  was  studied  using  three-dimensional  discrete  dislocation  dynamics 
(DDD).  A  systematic  interaction  model  between  dislocations  and  { 1 0 1 2}  tension  twin  boundaries  (TBs)  was  proposed  and  introduced  into  the 
DDD  framework.  In  addition,  a  nominal  grain  boundary  (GB)  model  based  on  experimental  results  was  also  introduced  to  mimic  the  GB’s  barrier 
effect.  The  current  simulations  show  that  tension  TBs  act  as  strong  obstacles  to  gliding  dislocations,  and  contribute  significantly  to  the  overall  hard¬ 
ening  response,  while  twin  growth  results  in  a  softening  effect.  Therefore,  the  Mg  concave  stress-strain  curve  can  be  explained  in  terms  of  the  com¬ 
petition  between  TB  induced  hardening  and  twin  growth  induced  softening.  At  low  strain  levels,  twin  growth  induced  softening  dominates  and  a 
decreasing  hardening  rate  is  observed  in  Stage-I.  In  Stage-II,  the  TB  induced  hardening  dominates,  which  leads  to  an  increasing  hardening  rate. 
©2015  Acta  Materialia  Inc.  Published  by  Elsevier  Ltd.  All  rights  reserved. 
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1.  Introduction 

With  the  increasing  demand  to  reduce  carbon  dioxide 
emission,  many  efforts  are  directed  toward  reducing  the 
structural  weight  of  fossil  fuel  powered  vehicles  to  increase 
their  fuel  efficiency.  As  the  lightest  structural  metal  (one- 
third  lighter  than  aluminum),  magnesium  (Mg)  and  its 
alloys  have  been  attracting  a  lot  of  attention  in  recent  years 
for  their  potential  use  in  automotive,  aerospace  and  defense 
applications  [1].  However,  at  present,  the  wide  use  of  Mg 
alloys  as  a  structural  material  is  still  challenging  due  to 
their  poor  room  temperature  formability.  Owing  to  its 
hexagonal  closed  packed  (hep)  lattice  structure  and  low 
crystal  symmetry,  Mg  cannot  maintain  sufficient  indepen¬ 
dent  slip  modes  to  accommodate  arbitrary  deformation  in 
poly  crystalline  materials  at  room  temperature. 
Consequently,  in  addition  to  dislocations,  twinning  plays 
an  important  role  in  their  plastic  deformation.  Therefore, 
it  is  vital  to  grasp  a  full  understanding  of  the  combined 
effects  of  dislocations,  twinning,  and  their  mutual 
interactions  to  address  their  role  on  the  overall  mechanical 
behavior  of  hep  materials  [2], 

Due  to  deformation  twinning,  the  deformation  of  Mg 
and  its  alloys  displays  characteristics  that  are  different  than 
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those  commonly  observed  in  face  centered  cubic  (fee)  and 
body  centered  cubic  (bcc)  crystals.  One  common  character¬ 
istic  is  that  the  true  stress-strain  curves  of  both  single  crys¬ 
tals  [3,4]  and  polycrystalline  [4,5]  Mg  display  a  concave 
shape.  In  addition,  the  stress-strain  curve  can  be  character¬ 
ized  by  three  distinct  stages.  In  Stage-I  a  decreasing  harden¬ 
ing  rate  is  observed,  followed  by  Stage-II  with  an  increasing 
hardening  rate,  and  finally  a  decreasing  hardening  rate  in 
Stage-Ill  [4,5],  Stage-I  is  generally  associated  with  twin- 
ning-mediated  plasticity  in  the  form  of  twin  nucleation 
and  growth,  leading  to  a  decreasing  hardening  rate  [6-8], 
On  the  other  hand,  the  progressive  decrease  in  the  work 
hardening  rate  in  Stage-Ill  is  a  result  of  post-saturation 
in  the  twin  volume  fraction  leading  to  predominant  disloca¬ 
tion-mediated  plasticity  [6,9].  In  Stage-II,  both  dislocation 
slip  and  twinning  are  active  and  interact;  however,  the 
underlying  mechanisms  associated  with  the  increasing 
hardening  rate  are  still  debated. 

Jiang  et  al.  [6]  examined  the  microstructure  evolution 
during  AZ31  deformation,  and  observed  that  in  Stage-II 
intersections  between  primary  and  secondary  { 1 0 1 2}  ten¬ 
sion  twins  result  in  significant  grain  refinement  and  subse¬ 
quently  an  increasing  hardening  rate.  On  the  contrary, 
Wang  et  al.  [10]  reported  that  the  barrier  effect  of  tension 
twin  boundaries  (TBs)  weakens  as  twins  coalesce  during 
Stage-II,  which  undermines  the  suggested  contribution 
from  grain  refinement.  In  addition,  they  suggested  that 
the  increasing  hardening  rate  can  be  attributed  to  texture 
induced  hardening.  This  is  a  result  of  the  tension  twin  hard 
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orientation,  since  its  c-axis  is  nearly  parallel  to  the  loading 
axis.  Knezevic  et  al.  [7]  showed  that  tension  twins  consume 
the  entire  grain  before  the  end  of  Stage-II,  and  all  the  afore¬ 
mentioned  twinning  mechanisms  cannot  explain  the  peak 
strain  hardening  rate  at  the  end  of  Stage-II.  Instead,  they 
suggested  that  this  peak  strain  hardening  rate  may  be 
attributed  to  the  nucleation  of  compression  twins  in  the 
tension  twinned  grains.  Barnett  et  al.  [5]  observed  the  for¬ 
mation  of  low  angle  boundaries  arising  from  the  glissile- 
to-sessile  transition  of  dislocations  induced  by  twinning 
shear  [11,12],  which  can  also  act  as  a  source  of  hardening. 
Finally,  de-twinning  as  a  result  of  strong  interactions 
between  dislocations  and  twins  was  also  suggested  to  con¬ 
tribute  to  the  observed  increasing  strain  hardening  rate  in 
Stage-II  [13].  These  conflicting  proposed  mechanisms  indi¬ 
cate  that  there  is  no  clear  consonance  yet  regarding  the 
mechanisms  responsible  for  the  increasing  strain  hardening 
rate  in  Stage-II. 

In  the  past  two  decades,  several  computational  methods 
have  been  developed  to  characterize  dislocation-mediated 
plasticity.  Of  those,  discrete  dislocation  dynamics  (DDD) 
is  one  of  the  most  efficient  methods  to  capture  disloca¬ 
tion-mediated  plasticity  at  the  micro  scale.  To  date,  several 
two-dimensional  [14,15]  and  three-dimensional  [16-19] 
DDD  methods  were  proposed  and  used  widely  to  study 
crystal  size  effects  (e.g.  [20-22]),  grain  size  effects  (e.g. 
[23-25])  and  intermittency  behavior  (e.g.  [26]).  While  most 
previous  DDD  studies  have  focused  on  fee  and  bcc  crystals, 
more  recently  several  studies  have  been  performed  on  hep 
crystals.  In  particular  DDD  simulations  of  hep  materials 
have  focused  on  dislocation  junction  formation  and 
strength  [27,28],  orientation  influence  on  the  grain  size 
effects  [25],  micro/nano-pillar  plasticity  [29,30],  elastic  ani¬ 
sotropy  [31],  and  the  role  of  plastic  anisotropy  on  micro¬ 
cracking  [32],  It  should  be  noted  that  in  these  studies  defor¬ 
mation  twinning  was  not  considered.  Nevertheless,  twin¬ 
ning  plays  an  important  and  sometimes  dominant  role  in 
the  mechanical  behavior  of  both  single  crystalline  and  poly¬ 
crystalline  Mg.  As  a  result,  DDD  simulations  not  account¬ 
ing  for  twinning  may  lead  to  inaccurate  predictions  of  the 
mechanical  behavior  of  hep  materials.  This  also  indicates 
the  importance  and  urgency  of  developing  a  dislocation- 
TB  interaction  model  within  the  framework  of  DDD  to 
account  for  deformation  twinning. 

As  aforementioned,  the  deformation  of  Mg  and  its 
alloys  is  of  great  interest  due  to  their  numerous  potential 
applications  in  industry.  However,  the  mechanisms  associ¬ 
ated  with  their  hardening  response,  mainly  in  Stage-II,  are 
not  well  characterized.  In  order  to  address  this,  we  propose 
a  systematic  interaction  model  between  dislocations  and 
{ 1 0 1 2}  TBs,  and  integrate  this  model  into  3D-DDD  sim¬ 
ulations.  A  simple  grain  boundary  (GB)  model,  which 
agrees  with  experimental  results,  is  also  implemented. 
Using  this  complex  DDD  framework,  the  deformation  of 
micro-twinned  polycrystalline  Mg  is  investigated. 
Specifically  these  simulations  identify  the  influences  of 
grain  orientation,  GBs,  TBs  and  twin  volume  fraction  on 
the  deformation  of  Mg  polycrystals. 

2.  Computational  methods 

In  the  following  simulations,  a  cubic  simulation  cell  is 
employed  to  model  a  representative  grain  in  a  bulk  mate¬ 
rial,  as  shown  schematically  in  Fig.  1(a),  where  periodic 


boundary  conditions  (PBCs)  are  imposed  along  all  three 
directions.  A  twin  lamella  having  a  thickness  /t  is  intro¬ 
duced  at  the  center  of  the  grain  with  size  /g.  Thus,  the  twin 
volume  fraction  is  given  by  ft  =  /t//g.  Here,  only  { 1 0 1 2} 
tension  twins  are  considered,  as  shown  by  the  relative  ori¬ 
entations  of  the  grain  and  twin  in  Fig.  1(a).  It  has  been 
shown  experimentally  that  the  number  of  twins  per  grain 
is  proportional  to  the  square  of  grain  size,  and  the  probabil¬ 
ity  of  having  more  than  one  twin  per  grain  decreases  con¬ 
siderably  for  grains  having  diameters  less  than  40  pm  [33], 
Thus,  given  the  grain  sizes  modeled  here,  the  current  pro¬ 
posed  simulation  cell  can  be  considered  as  a  representative 
twinned  grain  in  a  bulk  Mg  polycrystal  with  grain  size 
below  40  pm. 

In  Mg  polycrystals,  strain  hardening  can  be  attributed  to 
three  microstructural  features:  dislocations,  GBs  and  TBs. 
In  order  to  investigate  their  individual  contributions,  three 
additional  simplified  cubic  simulation  cells  are  considered, 
as  shown  in  Fig.  1.  In  one  case,  only  dislocations  are 
accounted  for  in  a  single  crystal  cubic  simulation  cell  with 
PBCs  along  all  three  directions.  In  another  case,  the  con¬ 
tribution  of  GBs  is  accounted  for  in  a  polycrystal  sim¬ 
ulation  cell  in  the  absence  of  TBs.  In  this  case,  the  six 
surfaces  of  the  simulation  cell  are  considered  as  representa¬ 
tive  GBs  with  PBCs  along  all  three  directions.  Finally,  the 
effect  of  TBs  is  accounted  for  in  a  twinned  crystal  sim¬ 
ulation  cell  in  the  absence  of  GBs.  Here,  a  twin  lamella  is 
inserted  in  the  middle  of  the  simulation  cell  and  PBCs  are 
employed  along  all  three  directions.  It  should  be  noted  that 
these  simulation  cells  are  idealized,  since  twinning  is  typi¬ 
cally  observed  experimentally  in  bulk  single  crystals  and 
polycrystals.  Nevertheless,  these  simulations  can  provide 
new  insights  into  the  individual  contributions  of  disloca¬ 
tions,  GBs  and  TBs  on  the  hardening  response  of  Mg. 

In  metals,  screw  dislocation  cross-slip  is  an  important 
mechanism  for  strain  hardening  [34]  and  creep  behavior 
[35].  Therefore,  a  general  cross-slip  mechanism  is  utilized 
in  the  current  simulations.  Cross-slip  of  a  screw  dislocation 
is  possible  when  more  than  one  slip  system  share  the  same 
Burgers  vector  [36],  In  particular  (a)  dislocations  can  cross¬ 
slip  between  basal,  prismatic  and  pyramidal  I  planes,  (c) 
dislocations  between  prismatic  slip  planes,  and  (c  +  a) 
dislocations  between  prismatic,  pyramidal  I  and  pyramidal 
II  planes.  Further  details  about  the  cross-slip  mechanism 
used  here  can  be  found  elsewhere  [36]. 

In  fully  tension-twinned  grains,  compressive  loading 
parallel  to  the  c-axis,  or  tensile  loading  perpendicular  to 
the  c-axis,  would  promote  compression  twinning  [7], 
These  compression  twins  appear  always  in  a  thin  band 
shape  with  low  volume  fraction  [7].  To  model  this  effect, 
in  a  small  number  of  simulations  we  also  introduce  a  nomi¬ 
nal  compressive  twin  with  a  volume  fraction  of  0.05  in  the 
center  of  the  fully  twinned  grain.  These  compression  TBs 
are  assumed  to  be  impenetrable  to  dislocations. 

All  current  simulations  were  conducted  using  the  3D- 
DDD  simulation  code,  ParaDiS  [16],  All  possible  disloca¬ 
tion  Burgers  vectors  and  slip  planes  in  hep  crystals,  which 
are  listed  in  Table  1,  were  accounted  for  in  these  sim¬ 
ulations.  It  should  be  noted  that  according  to  the  hep  lattice 
structure,  slip  of  (a)  dislocations  and  dislocations  with 
Burgers  vector  bt  (i.e.  twinning  dislocations)  on  the 
{101 2}  twinning  planes  are  also  possible.  Thus,  these 
dislocations  are  accounted  for  in  the  current  simulations 
as  well.  The  experimentally  measured  Peierls  stresses  for 
dislocations  on  the  basal  (0.52  MPa  [37]),  prismatic 
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(b)  Single  crystal 
(no  GBs  or  TBs) 


(c)  Poly  crystal 
(no  TBs) 


(d)  Twinned  crystal 
(no  GBs) 


Fig.  1.  Schematic  of  the  four  simulation  cells,  showing  twin  boundaries,  grain  boundaries,  and  orientations,  (a)  Twinned  polycrystal  with  both  GBs 
and  TBs;  (b)  single  crystal  without  GBs  or  TBs;  (c)  polycrystal  without  TBs;  and  (d)  twinned  crystal  without  GBs. 


(39.2  MPa  [38]),  and  pyramidal  planes  (105  MPa  [39])  were 
introduced  into  the  code.  These  values  are  also  in  good 
agreement  with  recent  molecular  dynamics  (MD)  cal¬ 
culations  [40-43].  Basic  Mg  parameters  used  in  the  current 
simulations  include:  shear  modulus,  G=  17  GPa;  Poisson 
ratio,  v  =  0.29;  magnitude  of  (a)  dislocation  Burgers 


vector,  b  =  0.325  nm;  axial  ratio,  da  =  1.6236;  and  mass 
density,  p  =  1738  kg  m~3. 

While  most  DDD  simulations  typically  adopt  a  linear 
dislocation  mobility  law,  recent  MD  simulations  suggest 
that  the  dislocation  mobility  in  Mg  is  nonlinear  [41].  At 
low  stress  levels  the  dislocation  velocity  was  shown  to 
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Table  1.  Slip  systems  introduced  in  the  current  DDD  framework.  The  underlined  four  symbols  indicate  the  slip  systems  of  initial  dislocations  in  the 
simulation  cell. 


Slip  planes 

Burgers  vectors 

(a)  (1120) 

(c)  (0001) 

(c  +  a)  (1123) 

b,  ;i[0 111]* 

Basal  -  {0001} 

X. 

Prismatic  -  {1010} 

X 

X 

Pyramidal  I  -  { 1 0 1 1 } 

>< 

X 

Pyramidal  II  -  {1  122} 

_x 

Twinning  plane  -  {1012} 

X 

X 

increase  linearly  with  the  applied  stress,  and  at  higher  stress 
levels  the  velocity  also  increases  linearly  but  with  a  much 
smaller  slope.  Thus,  in  the  current  simulations  we  employ 
an  MD  fitted  bi-linear  relationship.  The  drag  coefficients 
and  transition  velocities  for  the  stress-velocity  curve  are 
listed  in  Table  Al  in  Appendix  A.  The  upper  limit  to  the 
dislocation  velocity  was  set  to  be  the  shear  wave  velocity, 
yj G/p.  Furthermore,  as  a  first  approximation,  the  disloca¬ 
tion  mobility  and  Peierls  stress  on  the  twinning  plane  were 
set  equal  to  those  on  the  pyramidal  plane,  since  their  pre¬ 
cise  values  have  not  been  reported  in  literature  yet. 

Uniaxial  strain-controlled  loading  (compressive  or  ten¬ 
sile)  with  a  constant  strain  rate  of  e  =  5000  s_1  was  imposed 
in  the  current  simulations.  Due  to  the  configuration  sym¬ 
metry  along  the  x  and  z  axes  with  respect  to  the  twin 
lamella,  the  loading  axis  is  confined  within  the  yz  plane. 
As  a  result,  four  representative  loading  directions  denoted 
by  v,  z,  yz  and  yz,  are  considered.  Loadings  along  the  y 
and  z  directions  are  realized  by  imposing  a  constant  strain 
rate  ey  =  e  and  ez  =  e  in  the  y  and  z  directions,  respectively. 
The  yz  loading  is  along  a  direction  45°  from  either  the  y  or  z 
axis,  and  is  achieved  by  imposing  the  strain  rate  tensor: 


£ 


~-v  0  0 

e  0  (1  —  «)/2  (l  +  v)/2 

_0  (1  + 1>)/2  (l-p)/2 


(1) 


Similarly,  the  yz  loading  is  along  a  direction  45°  from  the 
y  axis  and  135°  from  the  z  axis,  and  is  achieved  by  imposing 
the  strain  rate  tensor: 


£ 


~-v  0  0 

e  0  (1  —  v)/2  —  ( 1  +v)/2 

_0  — (1  +  v)/2  (1  —  v)/2 


(2) 


Initially,  Frank-Read  (FR)  dislocation  sources  were 
randomly  distributed  within  the  simulation  cell  (in  both 
the  matrix  and  the  twinned  region).  The  initial  dislocation 
source  density  in  all  simulations  was  pSTC  =  5  x  1012  m  , 
and  the  FR  source  length  was  /src  =  0.26  pm.  While 
dislocation  slip  on  various  slip  planes  is  possible  in  hep 
crystals  as  listed  in  Table  1,  (a)  dislocations  on  basal, 
prismatic  and  pyramidal  I  planes  as  well  as  {c  +  a) 
dislocations  on  pyramidal  II  planes  play  a  predominant 
role  in  dislocation  mediated  plasticity  of  Mg  [44,45], 
Therefore,  in  the  current  simulations  all  dislocation 
sources  were  initially  assigned  to  these  four  slip  planes 
with  equal  numbers,  which  also  agree  with  previous  stud¬ 
ies  [25,46,47].  To  account  for  possible  statistical  variations 
due  to  the  random  distribution  of  the  dislocation  sources, 
each  simulation  was  repeated  three  times  with  different 
initial  distributions. 


2.1.  Dislocation-twin  boundary  interaction  model 

To  accurately  model  dislocation-TB  interactions  in 
DDD  we  first  derive  a  set  of  physically  based  rules  for  these 
interactions  from  experimental  observations  and/or  MD 
simulations.  Price  [48]  first  studied  these  interactions  in 
Zinc  and  showed  that  an  {a\)  dislocation,  which  intersects 
the  TB  in  a  screw  character,  can  transmit  through  the  TB 
without  any  remnants  on  the  TB.  On  the  other  hand, (<72} 
and  (fl3)  dislocations  in  the  matrix  can  only  be  transmitted 
mutually  through  the  TB  by  forming  a  (c)  dislocation  in  the 
twinned  crystal  while  leaving  two  twinning  dislocations  on 
the  TB.  Yoo  and  Wei  [49]  subsequently  developed  a  set  of 
dislocation-TB  interaction  scenarios  for  (a)  and  (c  +  a) 
dislocations  based  on  geometric  considerations  of  the  hep 
crystal  structure.  Electron  microscope  observations  have 
shown  good  agreement  with  the  Yoo  and  Wei  model  [49— 
52],  Furthermore,  Yoo  [53]  performed  a  detailed  analysis 
and  proposed  twenty  six  distinct  types  of  interactions 
between  three  dislocation  families  and  four  most  common 
twin  boundaries  in  hep  crystals.  In  addition  to  these  experi¬ 
mental  observations,  recent  MD  simulations  showed  that 
the  dislocation-TB  interactions  are  dependent  on  the  load¬ 
ing  direction.  In  particular,  it  was  shown  that  (ai)  disloca¬ 
tions  on  the  basal  plane  in  the  matrix  can  be  transmitted  to 
the  basal  plane  in  the  twin,  the  prismatic  plane  in  the  twin, 
or  the  prismatic  plane  in  the  matrix  depending  on  the  direc¬ 
tion  of  the  applied  load  [54,55].  In  addition,  no  residual  or 
twinning  dislocations  were  left  on  the  TB,  which  agrees  well 
with  experimental  observations.  However,  (a2)  and  (a3) 
dislocations  were  observed  to  introduce  a  twinning  disloca¬ 
tion  in  the  form  of  TB  steps  after  absorption  by  the  TB 
[56,57],  The  subsequent  motion  of  these  twinning  disloca¬ 
tions  leads  to  TB  migration,  which  agrees  with  the  twin 
growth  mechanism  by  twinning  dislocations  [56]. 

In  the  current  model,  two  rules  are  specified  to  identify 
the  reaction  outcomes  of  a  dislocation-TB  interaction,  in 
terms  of  twinning  dislocations  on  the  TB  as  well  as 
Burgers  vectors  and  slip  planes  of  the  transmitted  disloca¬ 
tions  in  the  twin.  In  particular,  these  rules  are  based  on 
(a)  geometric  constraints,  and  (b)  maximum  power 
dissipation. 

The  geometric  rule  predicts  the  resulting  dislocations 
and  twinning  dislocations.  Accordingly,  the  systematic 
interaction  scenarios  of  Yoo  [53]  between  different  disloca¬ 
tions  (i.e.  (a),  ( c )  and  (c  +  a))  and  the  {1012}  TB  were 
introduced  into  the  DDD  framework.  It  should  be  noted 
that  in  Yoo’s  model  some  peculiar  dislocation  reactions 
were  proposed.  In  particular,  (c  +  a)  and  (c)  incident 
dislocations  were  proposed  to  be  transmitted  as  disloca¬ 
tions  with  Burgers  vectors  [5  7  2  3]/6  and  [1 100],  respec¬ 
tively.  However,  these  high  lattice  index  dislocations 
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would  not  be  stable  because  of  their  high  Peierls  stresses 
and  core  energies,  and  should  be  decomposed  into  more 
common  dislocations  having  Burgers  vectors  with  low  lat¬ 
tice  indexes.  Our  MD  simulations  (not  presented  here) 
show  that  a  [1100]  dislocation  decomposes  into  two  (a) 
dislocations  such  that  [1 1 00]  — ►  [l2l0]/3  +  [2 1 1 0]/3  = 
(02)  —  (fli).  Furthermore,  simple  geometric  calculations 
show  that  [5723J/6  can  be  decomposed  into  1.5  (a)  and 
0.5  (c  +  a)  dislocations  (i.e.  [5723J/6 — ^  [1210J/2+ 

[2  1 1  3]/6  =  1.5(a)  +  0.5(c  +  a)).  Thus,  in  the  current  sim¬ 
ulations,  these  low  lattice  index  dislocations  are  defined 
instead. 

The  general  dislocation  decomposition  proposed  here  is 
listed  in  Table  A2  in  Appendix  A,  in  which  all  the  possible 
incident  dislocations  from  all  slip  planes  are  considered.  All 
the  resulting  dislocation  Burgers  vectors  and  slip  planes 
have  low  lattice  index.  The  resulting  slip  planes  were  chosen 
to  be  tilt  planes  with  respect  to  the  incident  plane.  Thus,  if  a 
tilt  plane  does  not  exist,  the  dislocation  is  not  allowed  to 
penetrate  (i.e.  twisted  boundaries  are  impenetrable).  This 
geometric  rule  is  analogous  to  that  used  for  dislocation- 
GB  interaction  in  DDD  simulations  of  fee  polycrystals 
[14,23,26]. 

While  for  most  incident  dislocations  this  geometric  rule 
is  sufficient  to  identify  the  transmitted  dislocation  planes, 
this  is  not  the  case  for  incident  (a  1)  dislocations,  which 
intersect  the  TB  in  a  screw  configuration.  In  this  case,  six 
potential  tilt  planes  in  the  twin  exist.  Serra  et.al  [54]  showed 
from  MD  simulations  that  depending  on  the  direction  of 
the  applied  load  a  matrix  (op)  dislocation  on  a  basal  plane 
can  be  transmitted  to  the  twin  basal  plane,  twin  prismatic 
plane,  or  reflected  back  into  the  matrix  on  a  prismatic 
plane.  In  these  simple  pure  shear  stress  simulations,  the 
resulting  slip  plane  with  the  largest  resolved  shear  stress 
was  typically  the  one  with  the  highest  dislocation  mobility. 
However,  in  the  current  DDD  simulations,  the  stress  field  is 
always  complex.  In  addition,  since  the  Peierls  stresses  and 
mobility  laws  are  strongly  dependent  on  the  slip  system 
under  consideration,  thus,  the  plane  with  largest  resolved 
shear  stress  does  not  necessarily  produce  the  highest 
dislocation  mobility.  To  account  for  all  these  factors,  the 
most  favorable  plane  is  chosen  to  be  the  one  that  would 
result  in  the  largest  power  dissipation  [16]: 

P='EFJVJ  (3) 

where  Fj  and  Vj  are  the  force  and  velocity  of  the  /th  trans¬ 
mitted  dislocation. 

It  should  be  noted  that  in  the  current  simulations,  plas¬ 
ticity  associated  with  TB  migration  or  twin  growth  will  not 
be  accounted  for  due  to  the  complexities  of  tracking  an 
evolving  boundary  within  the  framework  of  DDD. 

2.2.  Dislocation-grain  boundary  interaction  model 

For  the  representative  grain  model  with  PBCs  in  all 
three  directions,  shown  in  Fig.  1,  the  orientation  of  each 
grain  in  the  bulk  is  identical,  which  corresponds  to  a  strong 
basal  texture  in  polycrystalline  Mg.  However,  the  misori- 
entation  between  the  a-axes  in  neighboring  grains  is  not 
accounted  for  here.  Nevertheless,  in  the  current  sim¬ 
ulations,  the  GB  is  introduced  as  an  artificial  interface/bar¬ 
rier  to  dislocation  motion.  As  dislocations  approach  the 
GB  they  will  be  trapped.  As  a  result,  subsequent  disloca¬ 
tions  are  piled  up  at  the  GB,  leading  to  increasing  the  shear 


stress  on  the  leading  dislocation.  Once  this  shear  stress 
exceeds  the  GB  strength,  the  leading  dislocation  can  be 
transmitted  into  neighboring  grains.  This  approximation 
was  successfully  adopted  by  DDD  simulations  of  polycrys¬ 
talline  materials  [25,58],  The  GB  transmission  strength  is 
largely  governed  by  the  misorientation  angle,  0mis,  which 
can  be  approximated  as  rGB  =  2Gsin2(O.50mis)  [14,23,26], 
The  effect  of  the  misorientation  angle  on  the  stress-strain 
response  from  our  DDD  simulations  in  the  absence  of 
twinning  is  shown  in  Fig.  2.  In  these  simulations  the  grain 
size  was  /g  =  0.81  pm,  and  a  compressive  loading  along  the 
y-axis  was  imposed,  which  produces  a  predominant  basal 
slip  deformation.  An  experimental  stress-strain  curve  of  a 
0.8  pm  grained  Mg  alloy  under  a  similar  loading  condition 
[59]  is  also  shown.  It  is  clear  that  a  misorientation  angle  of 
0mis  =  15°  (i.e.  GB  strength  tGb  =  580  MPa)  agrees  well,  in 
terms  of  the  yield  strength,  with  this  experimental  result, 
and  thus  in  all  subsequent  simulations  this  GB  strength  is 
chosen. 


3.  Results  and  discussions 

3.1.  Hardening  response  of  single  crystals,  polycrystals  and 
twinned  crystals 

Here,  the  three  simplified  simulation  cells  shown  in 
Fig.  l(b-d)  are  employed  to  study  the  individual  con¬ 
tributions  of  dislocations  (i.e.  single  crystals),  GBs  (i.e. 
polycrystals  with  no  twins)  and  TBs  (i.e.  twinned  crystals) 
on  the  hardening  response  of  Mg.  All  the  simulation  cells 
have  an  edge-length  /g  =  1.3  pm,  and  in  the  simulations  of 
twinned  crystals  the  twin  thickness  is  /t  =  0.5/g.  Fig.  3  shows 
the  stress-strain  responses  under  different  loading  ori¬ 
entations  from  the  three  simulation  cells.  Fig.  4  shows  the 
contours  of  the  effective  plastic  strain  at  the  end  of  each 
simulation,  from  which  the  slip  traces  can  be  identified.  It 
is  observed  that  the  responses  of  single  crystals  loaded 
along  the  y  and  z  directions  (black  curves  in  Fig.  3)  are 
identical.  These  two  loading  directions  make  angles 


Fig.  2.  Engineering  stress-strain  responses  for  different  misorientation 
angles.  Experimental  result  from  [59]  is  also  shown  for  comparison. 
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Fig.  3.  Engineering  stress-strain  responses  of  single  crystals,  polycrys¬ 
tals  and  twinned  crystals  under  different  loading  directions.  Simulation 
cell  dimensions  are  lg=  1.3  pm,  and  /,  =  0.65  pm. 


46.85°  and  43.15°,  respectively,  with  respect  to  the  basal 
plane  leading  to  predominant  basal  slip,  as  shown  in 
Fig.  4(a),  since  the  Schmid  factor  on  these  planes  is  maxi¬ 
mum.  The  low  yield  stresses  for  these  two  loading  ori¬ 
entations  are  mainly  governed  by  the  low  Peierls  stress 
and  high  dislocation  mobility  on  the  basal  plane.  On  the 
other  hand,  loadings  along  the  yz  and  yz  directions  lead 
to  higher  yield  points.  The  yz  direction  is  almost  perpen¬ 
dicular  to  the  c-axis,  and  plasticity  is  mediated  by  prismatic 
(a)  dislocation  slip  (see  Fig.  4(b)).  Furthermore,  the  yz 
loading  is  1.85°  away  from  the  c-axis,  which  leads  to  pre¬ 
dominant  (c  +  a)  slip  on  pyramidal  II  planes  (see 
Fig.  4(c)).  As  expected,  the  yield  stress  for  loading  along 
the  yz  direction  is  the  highest  since  pyramidal  slip  is  the 
hardest  slip  mode.  It  is  also  observed  from  Fig.  3  that 
regardless  of  the  loading  directions,  all  single  crystals 
results  exhibit  weak  strain  hardening  response  in  the 
absence  of  GBs  or  TBs  up  to  1%  strain  for  this  selected 
small  periodic  simulation  cell. 

Fig.  3  also  shows  the  predicted  stress-strain  responses  of 
polycrystals  loaded  in  four  different  orientations  (red 
curves).  These  simulations  show  that  GBs  lead  to  a  signifi¬ 
cant  rise  in  the  crystal  strength  as  compared  to  single  crys¬ 
tal  simulations.  In  particular,  the  yield  strength  increases  by 
208%  for  y  and  z  loadings  (basal  slip),  by  70%  for  yz  load¬ 
ing  (prismatic  slip),  and  by  19%  for  yz  loading  (pyramidal 
slip),  which  are  directly  linked  to  the  GB’s  barrier  effect. 
From  the  contour  maps  in  Fig.  4(d-f)  the  deformation 
modes  are  observed  to  be  similar  to  those  from  single  crys¬ 
tal  simulations.  In  addition,  higher  strain  hardening  rates 
are  observed  in  all  simulations  with  GBs  (i.e.  the  mean  val¬ 
ues  from  three  realizations  with  different  initial  dislocation 
distributions  are  15.5  GPa  for  basal,  8.2  GPa  for  prismatic 
and  5.4  GPa  for  pyramidal  slips).  The  orientation  depen¬ 
dence  of  the  hardening  rate  is  a  result  of  the  anisotropy 
in  slip  modes  and  Peierls  stresses  in  Mg.  Since  basal  slip 
is  a  single  slip  orientation  (Fig.  4(d)),  many  dislocations  pile 
up  at  the  GBs,  which  lead  to  a  strong  back  stress  and  a  high 
strain  hardening  rate.  For  the  prismatic  and  pyramidal 
slip  orientations  (Fig.  4(e— f)),  plasticity  is  mediated  by 


multi-slip,  resulting  in  fewer  dislocation  pileups  and  lower 
strain  hardening  rates  at  the  same  strain  levels. 

Finally,  the  twinned  crystals  responses  are  shown  by 
blue  curves  in  Fig.  3.  A  twin  lamella  of  thickness  /t  =  /g/2 
is  introduced  at  the  center  of  the  simulation  cell.  For  these 
simulations,  only  two  unique  loading  directions  exist.  The  y 
and  z  loading  directions  share  the  same  stress-strain 
response  and  deformation  modes  in  the  matrix  and  twinned 
regions,  while  yz  and  yz  loading  directions  share  another 
one.  The  hardening  rates  for  the  y  and  yz  loading  directions 
are  20.4  GPa  and  26.6  GPa,  respectively,  for  deformation 
up  to  0.8%  strain.  These  rates  are  significantly  higher  than 
those  observed  from  the  polycrystal  simulations.  Fig.  4(g) 
shows  the  contours  of  the  effective  plastic  strain  for  loading 
along  the  y  direction.  Dislocation  slip  is  predominant  on 
basal  planes  in  both  the  matrix  and  the  twin.  Some  slip 
traces  are  observed  to  be  fully  blocked  in  the  twin  or 
matrix,  which  is  consistent  with  the  high  hardening  rate 
observed  in  Fig.  3.  On  the  other  hand,  some  traces  are  con¬ 
nected  across  the  TB,  forming  a  long  zigzag  trace  from 
basal  planes  in  the  matrix  to  basal  planes  in  the  twin. 
These  observations  are  in  qualitative  agreement  with 
experimental  observations  [49,50],  Furthermore,  Fig.  4(h) 
shows  a  different  deformation  response  during  loading 
along  the  yz  direction.  Plasticity  here  is  confined  in  the 
matrix  since  the  twin  is  in  a  much  harder  orientation  (pyr¬ 
amidal  slip)  as  compared  to  the  softer  orientation  of  the 
matrix  (prismatic  slip). 

In  order  to  identify  the  effects  of  TBs  on  the  hardening 
rate,  the  dislocation  configurations  from  the  twinned  crys¬ 
tal  simulations  in  both  loading  directions  are  shown  in 
Fig.  5.  For  simulations  with  loading  along  the  y  direction, 
it  is  observed  that  many  parallel  dislocations  pile  up  against 
the  TB  as  shown  in  Fig.  5(a)  by  straight  dislocations  (red  in 
the  matrix  and  green  in  the  twin).  In  this  case,  all  the  piled 
up  dislocations  in  the  matrix  and  twin  are  basal  {a2)  and 
(a3)  dislocations,  since  the  Schmid  factor  for  ( a t)  disloca¬ 
tion  is  zero.  According  to  the  dislocation-TB  interaction 
model  shown  in  Table  A2,  these  dislocations  should  decom¬ 
pose  as  follows 

{a2)  or  (a3)  — s-  0.5(e  +  a)  +  bt  (4) 

Thus,  when  the  first  dislocation  ( incident  dislocation  h ;) 
arrives  at  the  TB,  it  decomposes  into  one  half  (c  +  a) 
dislocation  (a  sessile  residual  dislocation  br  on  the  TB) 
and  a  twinning  dislocation  (At),  as  shown  in  Fig.  5(al). 
This  decomposition  is  similar  to  that  reported  by  MD  sim¬ 
ulations  [55-57]  and  predicted  experimentally  [49,50]. 
While  the  Schmid  factor  on  the  TB  is  zero,  the  twinning 
dislocation  is  observed  to  glide  away  on  the  TB,  as  shown 
in  Fig.  5(a2),  due  to  the  repulsive  force  between  the  residual 
dislocation  and  twinning  dislocation.  This  residual  disloca¬ 
tion  also  has  a  repulsive  interaction  with  subsequent  inci¬ 
dent  dislocations,  which  contributes  to  the  dislocations 
pileup  near  the  TB  (see  Fig.  5(a2)).  With  the  increasing 
applied  stress  and  number  of  pileup  dislocations,  the  shear 
stress  acting  on  the  new  leading  incident  dislocation  in  the 
pileup  increases  remarkably,  and  overcomes  the  repulsive 
interaction  with  the  residual  dislocation  on  the  TB.  Thus, 
the  new  incident  dislocation  would  reach  the  TB  and  then 
decompose  into  another  one  half  (c  +  a)  dislocation  and 
a  twinning  dislocation.  Therefore,  an  integrated  (c  +  a) 
dislocation  can  subsequently  be  emitted  into  the  twinned 
crystal  from  the  TB,  as  shown  in  Fig.  5(a3).  This  is 
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Fig.  4.  Contours  of  effective  plastic  strain  for:  single  crystals  with  loading  in  the  (a)  y,  ( b)  yz,  (c)  yz  directions;  polycrystals  with  loading  in  the  (d)  y, 
(e)  yz,  (f)  yz  directions;  and  twinned  crystals  with  loading  in  the  (g)  y,  (h)  yz  directions.  In  (g)  and  (h),  the  image  is  a  quad-duplicate  of  the  original 
contour. 


consistent  with  experimental  observations  [50,51].  As  fur¬ 
ther  dislocations  reach  the  TB,  the  same  decomposition 
processes  would  repeat.  Finally,  it  should  be  noted  that 
while  dislocation  transmission  across  the  TB  is  observed 
in  the  current  simulations,  it  occurs  only  at  high  local  stress 
levels  due  to  heavy  dislocation  pileups,  as  shown  in 
Fig.  5(a4).  This  is  one  possible  reason  for  the  strong  barrier 
effect  of  the  TB  based  on  the  current  proposed  dislocation- 
TB  interaction  model. 

Oil  the  other  hand,  during  loading  along  the  yz  direc¬ 
tion  two  dense  dislocation  networks  are  observed  on  the 


TBs  (see  Fig.  5(b)).  Fig.  5(bl)  shows  that  the  network  con¬ 
sists  of  two  intersecting  sets  of  parallel  straight  disloca¬ 
tions,  and  a  group  of  curved  dislocations.  As  stated 
previously,  during  loading  along  the  yz  direction,  plasticity 
is  accommodated  by  prismatic  slip  in  the  matrix. 
According  to  Eq.  (4),  the  straight  dislocations  are  half 
(c  +  a)  residual  dislocations,  resulting  from  the  incident 
prismatic  (a)  dislocations.  The  angle  between  the  intersect¬ 
ing  sets  is  45.7°,  which  is  the  angle  between  the  intersection 
lines  of  two  prismatic  planes  on  the  TB.  Furthermore,  the 
curved  dislocations  are  twinning  dislocations,  which  glide 
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Fig.  5.  Dislocation  configurations  in  twinned  crystals  loaded  in  the  (a)  y  and  (b)  y:  directions.  The  inserts  (al)  through  (a4)  and  (bl)  show 
dislocation-TB  interaction  details. 


on  the  TB,  since  the  Schmid  factor  is  nonzero  for  this 
orientation. 

As  shown  in  Fig.  4(g),  loading  along  the  y  direction  is  a 
single  slip  orientation  since  plasticity  is  accommodated  by 
basal  slip  only.  Thus,  plasticity  is  confined  to  distinct  planes 
and  many  dislocation  pileups  would  be  expected  at  the  TBs 
(see  Fig.  5(a)).  On  the  other  hand,  loading  along  the  yz 
direction  leads  to  multi-slip  deformation  on  prismatic 
planes.  Subsequently,  the  contour  of  effective  plastic  strain 
in  the  matrix  is  uniform,  as  seen  in  Fig.  4(h),  and  few 
dislocation  pileups  are  expected  at  the  TBs  (Fig.  5(b))  for 
the  same  strain  level.  It  should  be  noted  that  prismatic  slip 
during  yz  loading  induces  a  higher  hardening  rate  than 
basal  slip  during  y  loading,  which  seems  contradictory  to 
the  hardening  in  Mg  polycrystals  as  previously  discussed. 
In  fact,  this  is  a  result  of  the  hard  orientation  (pyramidal 
slip)  of  the  twin  with  respect  to  the  yz  loading  axis,  in  which 
almost  no  plastic  deformation  is  observed  (Fig.  4(h)). 
Furthermore,  in  the  yz  loading  the  Schmid  factor  on  the 
TB  plane  is  high,  leading  to  the  glide  of  many  twinning 


dislocations.  This  would  subsequently  contribute  greatly 
to  TB  motion  and  twin  growth  [56,57],  in  agreement  with 
experiments  oriented  for  compressive  loading  perpendicu¬ 
lar  to  the  c-axis  [6,8]. 

The  current  DDD  simulations  provide  direct  evidence 
of  the  three  popular  hardening  mechanisms  that  are  usu¬ 
ally  attributed  to  twinning  in  Mg  and  its  alloys:  (a) 
Hall-Petch  induced  hardening  due  to  grain  refinement 
by  the  evolving  twins  [9,60,61];  (b)  Basinski  mechanism: 
glissile-to-sessile  transformation  of  dislocations  due  to 
the  twinning  shear  [11,12];  and  (c)  texture  hardening 
induced  by  the  hard  orientation  in  the  twins  [10,12], 
Fig.  5(a)  and  (b)  shows  many  sessile  residual  disloca¬ 
tions  remain  on  the  TB,  which  agrees  well  with  the 
Basinski  mechanism.  The  TBs  also  act  as  strong  obsta¬ 
cles  to  dislocation  motion,  inducing  grain  refinement 
hardening.  In  addition,  the  absence  of  plasticity  in  the 
twinned  region  (see  Fig.  4(h))  clearly  demonstrates  the 
hardening  induced  by  the  lattice  orientation  of  the 
twins. 
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3.2.  Hardening  response  of  twinned  polycrystals 

To  identify  the  collective  contribution  of  dislocations, 
GBs  and  TBs  on  the  hardening  response  of  Mg,  in  this  sub¬ 
section  the  deformation  of  twinned  polycrystals  is  investi¬ 
gated.  In  the  following  simulations  the  grain  size  is  kept 
constant  at  /g  =  1.3  pm  and  the  twin  size,  /t,  is  varied  to 
address  the  effect  of  twin  volume  fraction.  Tension  twin¬ 
ning  has  been  reported  experimentally  to  be  facilitated  by 
compression  perpendicular  to  the  c-axis,  or  tension  parallel 
to  the  c-axis  [6,8].  These  correspond  to  the  yz  compressive 
loading  and  yz  tensile  loading,  respectively. 

The  engineering  stress-strain  responses  for  different  twin 
volume  fractions  (/[  =  /t//g)  under  yz  compressive  loading 
are  shown  in  Fig.  6(a).  The  case  of  f  =  0.0  corresponds 
to  a  polycrystal  prior  to  twinning,  while  /t=1.0  corre¬ 
sponds  to  a  fully  twinned  polycrystal.  It  is  clear  that  the 
yield  stress  increases  with  increasing  twin  volume  fraction, 
until  it  reaches  a  maximum  at  f  =  0.7.  The  inset  in  Fig.  6(a) 
shows  that  dislocation  activity  in  the  matrix  contributes 
mostly  to  the  overall  plastic  deformation  for  /t  <  0.7. 
Subsequently,  an  increasing  twin  volume  fraction  results 
in  reducing  the  matrix  volume,  which  in  turn  results  in 
increasing  the  yield  stress  of  the  crystal.  Furthermore,  the 
case  of  ft=  1.0  indicates  that  the  yield  stress  for  the  twin 
orientation  is  about  350  MPa.  Thus,  for  a  volume  fraction 
ft  >  0.5,  once  the  yield  stress  exceeds  350  MPa,  dislocation 
slip  in  the  twin  begins  to  contribute  to  the  overall  plastic 
deformation.  For  f  >  0.7,  the  overall  plastic  deformation 
is  mainly  accommodated  by  dislocation  slip  in  the  twin, 
since  fewer  dislocation  sources  would  be  present  in  the 
smaller  matrix  region.  Therefore,  the  yield  stress  decreases 
with  further  increase  in  the  twin  volume  fraction. 

In  a  fully  twinned  grain  (i.e./t  =  1.0),  the  c-axis  is  close 
to  the  yz  compressive  loading  direction,  which  is  an  ori¬ 
entation  that  would  promote  compression  twinning. 
Knezevic  et  al.  [7]  observed  compression  twins  and  double 
twins  in  such  oriented  grains.  These  compression  twins 
always  appear  in  a  thin  band  shape  having  a  low  volume 
fraction,  and  it  was  suggested  that  this  contributes  to  the 
peak  hardening  rate  in  Stage-II  [7],  To  model  this,  here, 
we  introduced  a  nominal  compressive  twin  with  a  volume 
fraction  of  0.05  in  the  center  of  the  fully  twinned  grain. 
In  the  current  simulations  these  compression  TBs  are 
assumed  to  be  impenetrable  to  dislocations.  To  keep  the 
notation  consistent,  this  case  is  denoted  by  /t=1.05. 
Fig.  6(a)  shows  that  the  volume  fraction  of f  =  1.05  results 
in  the  highest  yield  stress,  which  supports  the  potential  con¬ 
tribution  of  compression  twins  on  hardening  in  Mg  poly¬ 
crystals  after  tension  twin  saturation. 

The  stress-strain  responses  as  a  function  of  tension  twin 
volume  fraction  under  yz  tensile  loading  are  plotted  in 
Fig.  6(b).  The  yield  stress  displays  a  different  character  than 
that  of  yz  compressive  loading.  The  yield  stress  is  observed 
to  increase  with  increasing  twin  volume  fraction  up  to 
f  =  0.2.  Subsequently,  for  a  volume  fraction  0.2</t<  1.0, 
the  yield  stress  continuously  decreases.  This  can  be 
explained  as  follows.  Under  yz  tensile  loading,  the  matrix 
c-axis  is  close  to  the  loading  direction,  thus,  the  matrix  is 
oriented  for  a  hard  deformation  mode  (i.e.  pyramidal  slip 
orientation).  For  ft  ^  0.2,  the  plastic  deformation  is  mainly 
mediated  by  dislocations  activity  in  the  matrix.  Thus,  as  the 
twin  volume  fraction  increases  the  matrix  volume 
decreases,  leading  to  an  increasing  yield  stress.  For 


Fig.  6.  Engineering  stress-strain  responses  of  twinned  polycrystals 
having  different  twin  volume  fractions  under  (a)  yz  compressive  and  (b) 
yz  tensile  loadings. 

ft  >  0.2,  the  plastic  deformation  is  gradually  dominated 
by  dislocation  slip  in  the  twin,  which  is  in  a  soft  orientation 
(i.e.  prismatic  slip  orientation),  as  shown  by  the  inset  in 
Fig.  6(b).  As  a  result,  the  yield  stress  decreases  with  increas¬ 
ing  twin  volume  fraction. 

It  should  be  noted  that  the  initial  dislocation  densities  in 
all  simulations  of  different  twin  volume  fractions  shown  in 
Fig.  6  are  the  same.  However,  in  real  Mg  crystals,  the  twin 
volume  fraction  increases  with  increasing  strain,  and  it 
might  be  expected  that  subsequently  the  dislocation  density 
would  also  increase.  These  increases  in  the  dislocation  den¬ 
sity  are  not  accounted  for  in  the  current  study  due  to  com¬ 
putational  difficulties  associated  with  the  increasing 
number  of  degrees  of  freedom  with  increasing  dislocation 
density. 

Fig.  7  shows  the  hardening  rate  due  to  dislocation  slip, 
8S,  as  a  function  of  twin  volume  fraction  for  the  two  load¬ 
ing  directions.  These  rates  are  computed  from  a  linear  fit  of 
the  stress-strain  curves  in  Fig.  6  above  a  plastic  strain  of 
£p  =  —0.05%.  It  is  observed  that  the  dislocation  slip 
induced  hardening  rate  during  yz  loading  increases  with 


H.  Fan  el  al.  I  Acta  Materials  92  (2015)  126-139 


135 


increasing  twinning  volume  fraction  for  /t  ^  0.3,  then 
decreases  for  0.3  1.0,  as  shown  in  Fig.  7(a). 

Subsequently,  the  hardening  rate  increases  again  with  the 
introduction  of  a  compression  twin  with  0.05  volume  frac¬ 
tion  in  the  fully  twinned  grain  (i.e.  ft  =  1.05).  Oil  the  other 
hand,  the  dislocation  slip  induced  hardening  rate  during  yz 
tensile  loading  from  the  current  DDD  simulations  is  shown 
in  Fig.  7(b).  Similar  trends  to  the  yz  compressive  loading 
are  observed,  with  0S  increasing  for  ft  <  0.4,  then  decreas¬ 
ing  for  0.4  ^/t  ^  1.0.  Furthermore,  for  the  case  of  fully 
twinned  grains  with  compressive  twinning  at  a  volume  frac¬ 
tion  of  0.05  (i.e.  f  =  1.05),  the  dislocation  slip  induced 
hardening  rate  is  observed  to  subsequently  increase,  which 
agrees  well  with  experimental  predictions  [7], 

It  should  be  noted  that  only  dislocation-mediated  plas¬ 
ticity  is  accounted  for  in  the  current  DDD  simulations. 
Under  actual  experimental  conditions,  plastic  deformation 
is  mediated  by  both  dislocation  slip  and  twin  growth  (TG), 
with  twin  growth  playing  a  major  role  especially  at  low 
strain  levels.  Thus,  the  contribution  from  twinning  plastic¬ 
ity  must  be  taken  into  account.  Accordingly,  the  overall 
hardening  rate,  8,  can  be  decomposed  into  a  contribution 


Fig.  7.  Slip  induced  hardening  rate  0s  as  a  function  of  twin  volume 
fraction  under  (a)  yz  compressive  and  (b)  yz  tensile  loadings. 


induced  by  dislocation  slip  and  a  contribution  induced  by 
twin  growth  such  that: 


_  da  _  rf(£(e  ‘dip  Ci,,)) 
de  de 


Here,  0Tq  is  the  hardening  rate  induced  by  twin  growth 
only,  and  can  be  computed  from  the  twin  volume  fraction 
[8]  such  that 
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where  m  and  ytwin  are  the  average  Schmid  factor  of  the  twin 
variants  and  characteristic  twinning  shear,  respectively.  For 
{101 2}  tension  twins  in  Mg  and  its  alloys,  ytwin  =  0.13, 
while  m  =  0.437,  and  0.499  for  yz  compressive  loading 
and  for  yz  tensile  loading,  respectively  [8], 

For  comparison  with  experiments,  the  dislocation  slip 
induced  hardening  rates,  0s,  as  measured  from  AZ31  poly¬ 
crystals  experiments  under  yz  compressive,  and  yz  tensile 
loadings  [8]  are  also  shown  in  Fig.  7.  It  is  observed  that 
the  shape  of  the  dislocation  slip  induced  hardening  rate, 
8s,  from  DDD  simulations  is  in  qualitative  agreement  with 
the  experimental  measurements,  indicating  that  the  DDD 
simulations  are  capable  of  capturing  the  main  deformation 
mechanisms  in  Mg  polycrystals.  However,  the  experimental 
curves  show  higher  hardening  rates  than  the  current  DDD 
simulations.  This  difference  can  be  attributed  to  two  main 
reasons:  (1)  the  experiments  are  performed  on  AZ31  poly¬ 
crystals  having  a  high  density  of  second  phase  particles 
[29],  while  the  simulations  are  performed  on  pure  Mg  poly¬ 
crystals;  (2)  the  initial  dislocation  density  in  the  simulations 
is  always  fixed  at  pSIC  =  5  x  1012  m~2  for  all  modeled  twin 
volume  fractions,  while  in  the  experiments  the  dislocation 
density  is  expected  to  continuously  increase  with  increasing 
twin  volume  fraction  or  applied  strain.  These  differences 
between  the  simulation  conditions  and  the  experiments  sug¬ 
gest  that  the  dislocation-mediated  hardening  rate  as  mea¬ 
sured  experimentally  will  be  higher  than  those  predicted 
from  the  current  DDD  simulations. 

In  order  to  compute  the  twin  growth  induced  hardening 
rate  8TG  through  Eq.  (6),  the  twin  volume  fraction  reported 
experimentally  from  six  different  groups  [6-8,62-65]  for 
these  two  loading  paths  is  plotted  in  Fig.  8  as  a  function 
of  the  applied  strain.  It  is  important  to  note  that  the  twin 
volume  fraction  as  a  function  of  strain  has  been  reported 
in  literature  to  be  independent  of  grain  sizes  [66],  which 
agrees  well  with  Fig.  8  for  grain  sizes  10-107  pm. 
Through  this,  the  twin  growth  induced  hardening  rate, 
8 to,  can  be  numerically  computed  from  Eq.  (6)  based  on 
the  results  shown  in  Fig.  8. 

Finally,  Fig.  9  shows  the  dislocation  slip  induced  hard¬ 
ening  rate,  0s,  as  predicted  from  DDD,  the  twin  growth 
induced  hardening  rate,  0tg,  as  predicted  from  Fig.  8, 
and  the  overall  hardening  rate,  0  =  0s  +  0tg,  for  the  two 
simulated  loading  paths.  For  the  yz  compressive  loading, 
the  experimentally  measured  0TG  is  always  negative  and 
continuously  increases  with  increasing  twin  volume  frac¬ 
tion,  as  shown  in  Fig.  9(a).  Furthermore,  two  stages  can 
be  identified  from  the  overall  hardening  rate  0:  ( a)  decreas¬ 
ing  hardening  rate  for  twin  volume  fraction  ft  ^  0. 1;  and  (b) 
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Fig.  8.  Experimentally  reported  tension  twin  volume  fraction  in  AZ3 1 
with  different  grain  sizes  as  a  function  of  strain  in  two  different  loading 
paths.  The  best  polynomial  fitting  curves,  ft  =  — 42. Oe2  +  13. 6e  —  0.083 
for  yz  compressive  loading  (perpendicular  to  c-axis),  and 
/,  =  — 1650. 6e3  +  303.4e2  —  4.9e  +  0.00817  for  yz  tensile  loading  (par¬ 
allel  to  c-axis),  are  also  shown. 

increasing  hardening  rate  for  ft>  0.1.  These  two  stages 
mimic  the  experimentally  reported  Stage-I  and  Stage-II 
[6-8,64],  In  addition,  for  ft=  1.05  (i.e.  compression  twins 
in  fully  twinned  grains),  the  highest  hardening  rate  is 
achieved.  This  indicates  that  compression  twins  in  the  fully 
tension-twinned  grains  can  be  attributed  to  the  peak  hard¬ 
ening  rate  in  Stage  II,  which  agrees  well  with  experimental 
results  [7].  It  should  be  noted  that  the  total  hardening  rate 
is  mainly  negative,  because  the  dislocation  slip  induced  9S 
computed  in  current  DDD  simulations  is  lower  than  that 
measured  in  experiments,  as  clearly  shown  in  Fig.  7. 

As  previously  mentioned,  there  are  at  least  three  sources 
contributing  to  strain  hardening  in  Mg  polycrystals,  includ¬ 
ing  dislocation  forest  hardening,  GBs  induced  hardening, 
and  TBs  induced  hardening,  which  all  contribute  to  9S.  In 
Section  3.1,  the  dislocation  forest  hardening  effect  for  single 
crystals  was  observed  to  be  weak  (see  Fig.  3).  Furthermore, 
GBs  do  produce  considerable  hardening,  however,  this  con¬ 
tribution  does  not  change  with  increasing  twin  volume  frac¬ 
tion.  Therefore,  the  variations  in  9S  as  a  function  of  twin 
volume  fraction  are  largely  induced  by  the  TBs,  also  con¬ 
veying  TBs  produce  a  hardening  effect  on  the  overall 
response.  On  the  other  hand,  the  experimentally  measured 
twin  growth  induced  hardening  rate,  0TG,  is  always  nega¬ 
tive,  as  shown  in  Fig.  9(a),  indicating  twin  growth  produces 
a  softening  effect  in  the  presence  of  dislocation  slip.  As  a 
result,  the  two  stages  in  the  overall  hardening  rate,  9,  can 
be  explained  by  the  competition  between  TB  induced  hard¬ 
ening  and  twin  growth  induced  softening.  In  Stage-I,  the 
twin  growth  induced  softening  dominates,  so  a  decreasing 
overall  hardening  rate  is  seen.  In  Stage-II,  both  the  TB 
induced  hardening  and  twin  growth  induced  softening 
decline,  but  the  latter  declines  more  rapidly,  leading  to  a 
decreasing  overall  hardening  rate. 

On  the  other  hand.  Fig.  9(b)  shows  the  hardening  rates 
during  yz  tensile  loading,  in  which  identical  phenomenon  to 
the  yz  compressive  loading  is  observed.  The  overall  harden¬ 
ing  rate  6  displays  two  distinct  stages,  which  are  similar  to 


Fig.  9.  Hardening  rate  as  a  function  of  twin  volume  fraction  for  (a)  yz 
compressive  loading  and  (b)  yz  tensile  loading. 

yz  compressive  loading  and  experimental  observations 
[8,64],  Again,  the  case  of  fully  twinned  grains  with  compres¬ 
sive  twinning  at  a  volume  fraction  of  0.05  (i.e.  f  =  1.05) 
yields  the  peak  hardening  rate  in  Stage  II,  in  good  agree¬ 
ment  with  experimental  observations  [7],  These  two  stages 
can  also  be  interpreted  based  on  the  competition  between 
TB  induced  hardening  and  twin  growth  induced  softening. 


4.  Conclusions 

In  this  work,  a  systematic  interaction  model  between 
dislocations  and  { 1 0 T  2}  TBs  was  proposed.  In  this  model 
all  the  possible  incident  dislocations  were  taken  into 
account.  In  order  to  mimic  the  GB  strengthening,  a  nomi¬ 
nal  GB  model  agreeing  with  experimental  results  was  pro¬ 
posed.  Both  the  TB  and  GB  models  were  introduced  into 
the  DDD  framework  to  study  the  mechanical  behavior  of 
a  representative  twinned  Mg  grain.  It  was  shown  that 
TBs  are  strong  obstacles  to  dislocation  glide,  and  play  a 
dominant  role  in  the  hardening  response.  On  the  other 
hand,  the  plastic  deformation  induced  by  twin  growth  pro¬ 
duces  a  softening  effect.  Therefore,  a  competition  exists 
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between  TB  induced  hardening  and  twin  growth  induced 
softening.  This  competition  results  in  the  concave  stress- 
strain  curve.  At  low  strain  levels,  the  twin  growth  induced 
softening  dominates,  leading  to  a  decreasing  hardening  rate 
in  Stage-I.  In  Stage-II,  the  hardening  effect  dominates, 
which  leads  to  an  increasing  hardening  rate. 
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Table  Al.  The  bi-linear  mobility  law  parameters  used  in  the  current  simulations. 


Slip  plane 

Dislocation  character 

Drag  coefficient  below 
transition  velocity  ( 10  6  pa  s) 

Transition  velocity  (km/s) 

Drag  coefficient  above 
transition  velocity  (1CT6  pa  s) 

Basal  plane 

Edge 

4.68 

1.83 

81.25 

Screw 

13.60 

0.85 

1625.0 

Prismatic  plane 

Edge 

8.25 

1.61 

125.0 

Screw 

36.52 

0.80 

108.33 

Pyramidal  I  plane 

Edge 

116.07 

- 

- 

Screw 

- 

- 

- 

Table  A2.  Dislocation  decomposition  at  TB. 


Incident 
dislocation 
in  matrix 

Incident  plane 

Resulting 
dislocations  in 
twin 

Resulting 

planes 

Number  of 
twinning 
dislocations 

<ai>  = 

[2110J/3 

Ba  (0001) 

<ai>  =  [21 10]/3 

Depending  on 
the  power 

dissipation  rule. 

0 

Pr  (0110) 

Pyl  (0111) 

Pyl  (0111) 

Tp  (0112) 

Tp  (0112) 

<«2>  = 

[1210]/3 

Ba  (0001) 

0.5<c-ai>  = 

[2113]/6 

Pr  (0110) 

1 

Pr  (1010) 

Pyll  (2112) 

Pyl  (1011) 

Pyl  (lOll) 

Pyl  (1011) 

Pyl  (1101) 

Tp  (1012) 

Impenetrable 

Tp  (1012) 

Impenetrable 

<a3>  = 

[Il20]/3 

Ba  (0001) 

-0.5<c+ai>  = 

[2113J/6 

Pr  (0110) 

-1 

Pr  (1100) 

Pyll  (2112) 

Pyl  (1101) 

Pyl  (1011) 

Pyl  (1101) 

Pyl  (1101) 

Tp  (1102) 

Impenetrable 

Tp  (1102) 

Impenetrable 

<c+a\>  - 

[2113]/3 

Pr  (0110) 

-2<«3>  = 

2[1120]/3 

Ba  (0001) 

-2 

Pyl  (lOTI) 

Pyl  (1101) 

Pyl  (1101) 

Pyl  (lTOl) 

Pyll  (2112) 

Pr  (1100) 

<c+a2>  = 

[1213]/3 

Pr  (1010) 

0.5<c-a3>  + 

1.5<«2>  = 

[1123]/6  + 

[1210J/2 

Pyll  (1122) 

Pr  (1010) 

-1 

Pyl  (0111) 

Pyl  (0111) 
Ba  (0001) 

Pyl  (1101) 

Pyl  (lOll) 
Pyl  (1011) 

Pyll  (1212) 

Pr  (1100) 

Tp  (1012) 
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<c+a3>  = 

[Il23]/3 

Pr  (1100) 

-0.5<c+a3>  + 

1.5<a2>  = 

[1123J/6  + 

[1210J/2 

Pr  (1100) 

Tp  (1012) 

-3 

Pyl  (0111) 

Pyl  (0111) 
Ba  (0001) 

Pyl  (1011) 

Pyl  (1011) 
Pyl  (1011) 

Pyll  (1122) 

Pyll  (1122) 
Impenetrable 

<c — fl  i  >  = 

[2113]/3 

Pr  (0110) 

2<a2>  =  2[1210]/3 

Ba  (0001) 

-2 

Pyl  (1011) 

Pyl  (1011) 

Pyl  (1101) 

Pyl  (1011) 

Pyll  (2112) 

Pr  (1010) 

<c-ai>  - 

[1213]/3 

Pr  (1010) 

-0.5<c— <22>  - 

1 .5<«3>  = 

[12l3]/6  + 

[1120J/2 

Pr  (1010) 

Tp  (1102) 

-3 

Pyl  (01H) 

Pyl  (0111) 
Ba  (0001) 

Pyl  (1101) 

Pyl  (1101) 
Pyl  (1101) 

Pyll  (1212) 

Pyll  (1212) 
Impenetrable 

<c-ai>  = 

[1123]/3 

Pr  (1100) 

0.5<c+a2>  - 

1.5<fl3>  = 

[1213]/6  + 

[1120]/2 

Pyll  (1212) 

Pr  (1100) 

-1 

Pyl  (0111) 

Pyl  (0111) 
Ba  (0001) 

Pyl  (10U) 

Pr  (1101) 

Pr  (1101) 

Pyll  (1122) 

Pyl  (1010) 
Tp  (1102) 

<c>  = 

[0001] 

Pr  (1100) 

<fl2>-<£Z3>  = 

[12l0]/3  + 

[1120J/3 

Pr  (1100) 

Tp  (1012) 

-2 

Pr  (0110) 

Ba  (0001) 

Ba  (0001) 

Pr  (1010) 

Tp  (1102) 

Pr  (1010) 

Ba  =  basal  plane,  Pr  =  prismatic  plane,  Pyl  =  pyramidal  1  plane, 
Pyll  =  pyramidal  II  plane,  Tp  =  twinning  plane. 


References 

[1]  A. A.  Luo,  J.  Magnesium  Alloy  1  (2013)  2. 

[2]  J.W.  Christian,  S.  Mahajan,  Prog.  Mater.  Sci.  39  (1995)  1. 

[3]  K.E.  Prasad,  K.  Raiesh,  U.  Ramamurty,  Acta  Mater.  65 
(2014)  316. 

[4]  E.W.  Kelley,  W.F.  Elosford,  Technical  Report,  1967. 

[5]  M.R.  Barnett,  Z.  Keshavarz,  A.G.  Beer,  D.  Atwell,  Acta 
Mater.  52  (2004)  5093. 

[6]  L.  Jiang,  J.J.  Jonas,  A.A.  Luo,  A.K.  Sachdev,  S.  Godet, 
Mater.  Sci.  Eng.  A  445^146  (2007)  302. 

[7]  M.  Knezevic,  A.  Levinson,  R.  JJarris,  R.K.  Mishra,  R.D. 
Doherty,  S.R.  Kalidindi,  Acta  Mater.  58  (2010)  6230. 

[8]  S.-G.  Hong,  S.H.  Park,  C.S.  Lee,  Acta  Mater.  58  (2010)  5873. 

[9]  A.A.  Salem,  S.R.  Kalidindi,  R.D.  Doherty,  Acta  Mater.  51 
(2003)  4225. 

[10]  B.  Wang,  R.  Xin,  G.  Huang,  Q.  Liu,  Mater.  Sci.  Eng.  A  534 
(2012)  588. 

[11]  Z.S.  Basinski,  M.S.  Szczerba,  M.  Niewczas,  J.D.  Embury,  S.J. 
Basinski,  Rev.  Metall.  94  (1997)  1037. 

[12]  A.A.  Salem,  S.R.  Kalidindi,  R.D.  Doherty,  S.L.  Semiatin, 
Metall.  Mater.  Trans.  A  37  (2006)  259. 

[13]  D.  Sarker,  D.L.  Chen,  Scr.  Mater.  67  (2012)  165. 

[14]  H.  Fan,  Z.  Li,  M.  Huang,  X.  Zhang,  Int.  J.  Solids  Struct.  48 
(2011)  1754. 

[15]  E.  Van  der  Giessen,  A.  Needleman,  Model.  Simul.  Mater.  Sci. 
Eng.  3  (1995)  689. 

[16]  A.  Arsenlis,  W.  Cai,  M.  Tang,  M.  Rhee,  T.  Oppelstrup,  G. 
Hommes,  T.G.  Pierce,  V.V.  Bulatov,  Model.  Simul.  Mater. 
Sci.  Eng.  15  (2007)  553. 


[17]  H.M.  Zbib,  M.  Rhee,  J.P.  Hirth,  Int.  J.  Mech.  Sci.  40  (1998) 
113. 

[18]  L.P.  Kubin,  G.  Canova,  M.  Condat,  B.  Devincre,  V.  Pontikis, 
Y.  Brechet,  Solid  State  Phenom.  23  and  24  (1992)  455. 

[19]  J.A.  El-Awady,  S.  Bulent  Biner,  N.M.  Ghoniem,  J.  Mech. 
Phys.  Solids  56  (2008)  2019. 

[20]  H.  Fan,  Z.  Li,  M.  Huang,  Scr.  Mater.  67  (2012)  225. 

[21]  C.  Ouyang,  Z.  Li,  M.  Huang,  H.  Fan,  Int.  J.  Solids  Struct.  47 
(2010)  3103. 

[22]  J.A.  El-Awady,  Nat.  Commun.  6  (2015). 

[23]  Z.  Li,  C.  Hou,  M.  Huang,  C.  Ouyang,  Comput.  Mater.  Sci.  46 
(2009)  1124. 

[24]  M.  Huang,  Z.  Li,  J.  Tong,  Int.  J.  Plast.  61  (2014)  112. 

[25]  H.  Fan,  S.  Aubry,  A.  Arsenlis,  J.A.  El-Awady,  Scr.  Mater.  97 
(2015)  25. 

[26]  H.  Fan,  Z.  Li,  M.  Huang,  Scr.  Mater.  66  (2012)  813. 

[27]  C.C.  Wu,  P.W.  Chung,  S.  Aubry,  L.B.  Munday,  A.  Arsenlis, 
Acta  Mater.  61  (2013)  3422. 

[28]  G.  Monnet,  B.  Devincre,  L.P.  Kubin,  Acta  Mater.  52  (2004) 
4317. 

[29]  Z.H.  Aitken,  H.  Fan,  J.A.  El-Awady,  J.R.  Greer,  J.  Mech. 
Phys.  Solids  76  (2015)  208. 

[30]  G.S.  Kim,  Ph.D.  thesis,  University  of  Grenoble,  2011. 

[31]  L.  Capolungo,  I.J.  Beyerlein,  Z.Q.  Wang,  Model.  Simul. 
Mater.  Sci.  Eng.  18  (2010)  085002. 

[32]  S.  Olarnrithinun,  S.S.  Chakravarthy,  W.A.  Curtin,  J.  Mech. 
Phys.  Solids  61  (2013)  1391. 

[33]  M.R.  Barnett,  Scr.  Mater.  59  (2008)  696. 

[34]  A.M.  Hussein,  S.L  Rao,  M.D.  Uchic,  D.M.  Dimiduk,  J.A.  El- 
Awady,  Acta  Mater.  85  (2015)  180. 


H.  Fan  el  al.  I  Acta  Materialia  92  (2015)  126-139 


139 


[35]  J.  Chevy,  F.  Louchet,  P.  Duval,  M.  Fivel,  Philos.  Mag.  Lett. 
92  (2012)  262. 

[36]  S.  Aubry,  V.V.  Bulatov,  M.  Rhee,  G.  Hommes,  A.  Arsenlis, 
Model.  Simul.  Mater.  Sci.  Eng.  (submitted  for  publication 
(2015)). 

[37]  H.  Conrad,  W.D.  Robertson,  Trans.  Met.  Soc.  AIME  209 
(1957)  503. 

[38]  R.E.  Reed-Hill,  W.D.  Robertson,  Acta  Metall.  5  (1957)  717. 

[39]  A.  Staroselsky,  L.  Anand,  Int.  J.  Plast.  19  (2003)  1843. 

[40]  Y.  Tang,  J.A.  El-Awady,  Mater.  Sci.  Eng.  A  618  (2014)  424. 

[41]  S.  Groh,  E.B.  Marin,  M.F.  Horstemeyer,  D.J.  Bammann, 
Model.  Simul.  Mater.  Sci.  Eng.  17  (2009)  075009. 

[42]  J.A.  Yasi,  T.  Nogaret,  D.R.  Trinkle,  Y.  Qi,  L.G.  Hector  Jr., 
W.A.  Curtin,  Model.  Simul.  Mater.  Sci.  Eng.  17  (2009) 
055012. 

[43]  H.  Fan,  J.A.  El-Awady,  Q.  Wang,  J.  Nucl.  Mater.  458  (2015) 
176. 

[44]  Y.  Tang,  J.A.  El-Awady,  Acta  Mater.  71  (2014)  319. 

[45]  S.R.  Agnew,  O.  Duygulu,  Int.  J.  Plast.  21  (2005)  1161. 

[46]  J.  Zhang,  S.P.  Joshi,  J.  Mech.  Phys.  Solids  60  (2012)  945. 

[47]  J.  Zhang,  K.T.  Ramesh,  S.P.  Joshi,  Model.  Simul.  Mater.  Sci. 
Eng.  22  (2014)  075003. 

[48]  P.B.  Price,  Proc.  R.  Soc.  Lond.  A  Mater.  260  (1961)  251. 

[49]  M.H.  Yoo,  C-T.  Wei,  Philos.  Mag.  14  (1966)  573. 

[50]  D.I.  Tomsett,  M.  Bevis,  Philos.  Mag.  19  (1969)  129. 


[51]  J.I.  Dickson,  C.  Robin,  Mater.  Sci.  Eng.  11  (1973)  299. 

[52]  S.  Lay,  G.  Nouet,  Philos.  Mag.  A  70  (1994)  1027. 

[53]  M.H.  Yoo,  Trans.  Met.  Soc.  AIME  245  (1969)  1841. 

[54]  A.  Serra,  D.J.  Bacon,  Acta  Metall.  43  (1995)  4465. 

[55]  M.  Yuasa,  K.  Masunaga,  M.  Mabuchi,  Y.  Chino,  Philos. 
Mag.  94  (2013)  285. 

[56]  A.  Serra,  D.J.  Bacon,  Philos.  Mag.  A  73  (1996)  333. 

[57]  J.  Wang,  I.J.  Beyerlein,  C.N.  Tome,  Int.  J.  Plast.  56  (2014) 
156. 

[58]  C.  Zhou,  R.  LeSar,  Int.  J.  Plast.  30-31  (2011)  185. 

[59]  B.  Li,  S.  Joshi,  K.  Azevedo,  E.  Ma,  K.T.  Ramesh,  R.B. 
Figueiredo,  T.G.  Langdon,  Mater.  Sci.  Eng.  A  517  (2009)  24. 

[60]  S.  Asgari,  E.  El-Danaf,  S.  Kalidindi,  R.  Doherty,  Metall. 
Mater.  Trans.  A  28  (1997)  1781. 

[61]  Y.  Xin,  M.  Wang,  Z.  Zeng,  M.  Nie,  Q.  Liu,  Scr.  Mater.  66 
(2012)  25. 

[62]  H.J.  Yang,  S.M.  Yin,  C.X.  Huang,  Z.F.  Zhang,  S.D.  Wu, 
S.X.  Li,  Y.D.  Liu,  Adv.  Eng.  Mater.  10  (2008)  955. 

[63]  C.  Lou,  X.  Zhang,  G.  Duan,  J.  Tu,  Q.  Liu,  J.  Mater.  Sci. 
Technol.  30  (2014)  41. 

[64]  J.H.  Park,  H.L.  Kim,  J.E.  Jung,  Y.W.  Chang,  Met.  Mater. 
Int.  19  (2013)  389. 

[65]  A.J.  Levinson,  Ph.D.  thesis,  2012. 

[66]  A.  Ghaderi,  M.R.  Barnett,  Acta  Mater.  59  (2011)  7824. 


